Hands-on Exercise 11: Calibrating Spatial Interaction Models with R

1.0 Introduction

1.1 Getting Started

In this hands-on exercise, we will learn how to calibrate SIM to determine factors affecting the public bus passenger flows during the morning peak in Singapore.

1.2 Overview

Spatial Interaction Models (SIMs) are mathematical models for estimating flows between spatial entities developed by Alan Wilson in the late 1960s and early 1970, with considerable uptake and refinement for transport modelling.

There are four main types of traditional SIMs:

• Unconstrained

• Production-constrained

• Attraction-constrained

• Doubly-constrained

1.3 Installing and loading R packages

pacman::p_load(tmap, sf, sp,
               performance, reshape2,
               ggpubr, tidyverse,
               DT, stplanr)

2.0 Data Acquisition

We will be using 2 datasets in this exercise:

• Weekday morning peak passenger flows at planning subzone level (i.e. od_data.rds)

• URA Master Plan 2019 Planning Subzone boundary in simple feature tibble data frame format (i.e. mpsz.rds)

3.0 Geospatial and Aspatial Data Handling

We will be using the st_read() from sf package to import the data into RStudio.

3.1 Importing Geospatial Data

mpsz <- st_read("data/geospatial",
                layer = "MPSZ-2019")

Reading layer `MPSZ-2019' from data source 
  `C:\kt526\IS415-GAA\Hands-on_Ex\Hands-on_Ex11\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

mpsz_sp <- as(mpsz, "Spatial")
mpsz_sp

class       : SpatialPolygonsDataFrame 
features    : 332 
extent      : 103.6057, 104.0885, 1.158699, 1.470775  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs 
variables   : 6
names       : SUBZONE_N, SUBZONE_C, PLN_AREA_N, PLN_AREA_C,       REGION_N, REGION_C 
min values  : ADMIRALTY,    AMSZ01, ANG MO KIO,         AM, CENTRAL REGION,       CR 
max values  :    YUNNAN,    YSSZ09,     YISHUN,         YS,    WEST REGION,       WR 

Computing the distance matrix

dist <- spDists(mpsz_sp, 
                longlat = FALSE)

sz_names <- mpsz$SUBZONE_C

colnames(dist) <- paste0(sz_names)
rownames(dist) <- paste0(sz_names)

3.2 Importing Aspatial Data

population <- read.csv("data/aspatial/pop.csv")

4.0 Preparing Flow Data

4.1 Importing the Origin Destination data

odbus <- read_csv("data/aspatial/origin_destination_bus_202210.csv")

odbus$ORIGIN_PT_CODE <- as.factor(odbus$ORIGIN_PT_CODE)
odbus$DESTINATION_PT_CODE <- as.factor(odbus$DESTINATION_PT_CODE) 

4.1.1 Extracting the study data

odbus6_9 <- odbus %>%
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 6 &
           TIME_PER_HOUR <= 9) %>%
  group_by(ORIGIN_PT_CODE,
           DESTINATION_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))

datatable(odbus6_9)

Writing to rds

write_rds(odbus6_9, "data/rds/odbus6_9.rds")

4.2 Importing bus stop data

busstop <- st_read(dsn = "data/geospatial",
                   layer = "BusStop") %>%
  st_transform(crs = 3414)

Reading layer `BusStop' from data source 
  `C:\kt526\IS415-GAA\Hands-on_Ex\Hands-on_Ex11\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21

4.3 Importing MPSZ-2019

mpsz <- st_read(dsn = "data/geospatial",
                   layer = "MPSZ-2019") %>%
  st_transform(crs = 3414)

Reading layer `MPSZ-2019' from data source 
  `C:\kt526\IS415-GAA\Hands-on_Ex\Hands-on_Ex11\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

Write data to rds

mpsz <- write_rds(mpsz, "data/rds/mpsz.rds")

4.4 Combing bus stop and mpsz

busstop_mpsz <- st_intersection(busstop, mpsz) %>%
  select(BUS_STOP_N, SUBZONE_C) %>%
  st_drop_geometry()

datatable(busstop_mpsz)

Write data to rds

write_rds(busstop_mpsz, "data/rds/busstop_mpsz.rds")  

od_data <- left_join(odbus6_9 , busstop_mpsz,
            by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C,
         DESTIN_BS = DESTINATION_PT_CODE)

To check for duplicates

duplicate <- od_data %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

od_data <- unique(od_data)

od_data <- left_join(od_data , busstop_mpsz,
            by = c("DESTIN_BS" = "BUS_STOP_N")) 

duplicate <- od_data %>%
  group_by_all() %>%
  filter(n()>1) %>%
  ungroup()

od_data <- unique(od_data)